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Feynman and Hibbs were the first to variational^ determine an effective potential 
whose associated classical canonical ensemble approximates the exact quantum partition 
function. We examine the existence of a map between the local potential and an 
effective classical potential which matches the exact quantum equilibrium density and 
partition function. The usefulness of such a mapping rests in its ability to readily 
improve Born-Oppenheimer potentials for use with classical sampling. We show that 
such a map is unique and must exist. To explore the feasibility of using this result to 
improve classical molecular mechanics, we numerically produce a map from a library of 
randomly generated one-dimensional potential/effective potential pairs then evaluate its 
performance on independent test problems. We also apply the map to simulate liquid 
para-hydrogen, finding that the resulting radial pair distribution functions agree well with 
path integral Monte Carlo simulations. The surprising accessibility and transferability of 
the technique suggest a quantitative route to adapting Born-Oppenheimer potentials, with 
a motivation similar in spirit to the powerful ideas and approximations of density functional 
theory. 

Keywords: effective potentials, path integral molecular dynamics, nuclear quantum propagation, liquid hydrogen, 
density functional theory 



INTRODUCTION 

The energy and mass scales of chemical motion lie in a regime 
between quantum and classical mechanics but for reasons of com- 
putational complexity, molecular modeling (MM) is largely per- 
formed according to Newton's laws. When classical Hamiltonians 
are chosen to reproduce properties of real material, classical 
MM is an efficient compromise. An increasing amount of MM 
uses highly accurate Born-Oppenheimer (BO) potential energy 
surfaces, which allow one to study complex bond rearrange- 
ments where experiment cannot motivate a potential (Car and 
Parrinello, 1985; Wang et al, 2013). The BO surface is incompat- 
ible with classical statistical mechanics in the sense that we would 
not expect a classical simulation on the BO surface to repro- 
duce properties of the real material, except in the limit of infinite 
temperature. 

Many approaches already exist to bridge this gap and 
study quantum equilibrium properties and dynamics: path inte- 
gral Monte Carlo (PIMC), ring polymer molecular dynamics 
(RPMD), centroid molecular dynamics (CMD), variational path- 
integral approximations, discretized path-integral approxima- 
tions, semi-classical approximations, thermal Gaussian molecular 
dynamics and colored-noise thermostats (Whitlock et al, 1979; 
Chandler and Wolynes, 1981; Jang et al., 2001; Nakayama and 
Makri, 2003; Poulsen et al, 2003; Liu and Miller, 2006; Paesani et 
al, 2006; Fanourgakis et al, 2009; Liu et al., 2009; Ceriotti et al, 
201 1; Georgescu et al., 201 1; Perez and Tuckerman, 201 1 ). Most of 
the these methods involve computational overhead significantly 
beyond classical mechanics and as they approach exactness their 
cost rapidly increases. 



An alternative philosophy is suggested by density functional 
theory (DFT) (Hohenberg and Kohn, 1964; Mermin, 1965; Sham 
and Kohn, 1966; Runge and Gross, 1984; Yuen-Zhou et al, 2010; 
Tempel and Aspuru-Guzik, 2011). Following this line of reason- 
ing, three questions arise. Can an equilibrium quantum density 
be obtained from purely classical mechanics and an effective 
Hamiltonian? Is the effective Hamiltonian uniquely determined 
by the physical potential? Can the particle density and free energy 
be given by such a fictitious system? To all these questions, the 
answer "yes" is implied by the usual recipe for classical force-fields 
that fit experimental data. The idea of using an effective classical 
Hamiltonian to incorperate nuclear quantum propagation effects 
is not novel. For the first time, we prove the uniqueness and exis- 
tence of a map yielding a classical effective potential given the 
physical potential. We also make a contribution to this field by 
demonstrating that the aforementioned mapping can be reversed 
numerically and approximated analytically. 

The bargain of our proposed effective classical potential is 
similar to that posed by DFT. One sacrifices access to rigorous 
momentum based-observables and abandons the route to sys- 
tematic improvement. In exchange, the two properties which are 
physically guaranteed, the equilibrium particle density and the 
partition function, are obtained at a cost equivalent to classical 
sampling but with improved accuracy. As a practical tool, the 
map is an easy way to transform BO-based force fields into a 
form which is well-suited for classical sampling. Perhaps the most 
promising aspect of this mapping would be its scalability which 
could potentially extend the ability to treat quantum propaga- 
tion effects to all systems that can be sampled classically. It is 
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even possible that the fictitious trajectories of particles moving on 
such a potential would, like Kohn-Sham orbitals, have somewhat 
improved physicality over their classical counterparts, although 
we will not examine that possibility here. 

First, we show the uniqueness of an equilibrium effective 
potential that gives the exact equilibrium quantum density via 
classical sampling. Next, we demonstrate that the equilibrium 
effective potential may be approximated by a linear operator act- 
ing on the true potential. Finally, we numerically approximate the 
map in a rudimentary way, and obtain surprisingly good results 
and transferability for both one dimensional potentials and a 
model of liquid para-hydrogen. 

1. EQUILIBRIUM EFFECTIVE POTENTIAL 

In their seminal work on path integral quantum mechanics, 
Feynman and Hibbs introduced the concept of an effective classi- 
cal potential that allows for the calculation of quantum partition 
functions in a seemingly classical fashion (Feynman and Hibbs, 
1965). In Appendix A, we discuss a connection with the large 
and fruitful body of research that focuses on the centroid effec- 
tive potential and density which should not be confused with the 
equilibrium effective potential that we now define (Giachetti and 
Tognetti, 1985; Feynman and Kleinert, 1986; Voth, 1991; Cao, 
1993, 1994a,b,c; Cuccoli et al, 1995; Cao and Martyna, 1996; 
Martyna, 1996; Pavese and Voth, 1996; Roy et al., 1999; Krajewski 
and Muser, 2001; Blinov and Roy, 2004; Hone et al, 2005; Roy, 
2005; Mielke et al, 2013). We start by considering the equilibrium 
density matrix, 



(i) 



where TL is the system Hamiltonian, (3 is the inverse temperature, 
and Z is the partition function. Feynman showed us that we could 
connect this expression to the path integral representation of the 
quantum propagator 1 , 



-AMx)) 



Po(q a ,qb) = -Ty I Vr(x)e 

Z Jr(0) = q a 

where the Wick- rotated (t — > —ix) action functional is. 

N 

J2 ^n(t) 2 + V(r(x)) 



(2) 



1 fV' 
A[r(x)] = -J^ dx 



(3) 



By integrating over only closed paths at each coordinate we obtain 
the scalar equilibrium density, 



I I rr($h) = q 

Tlo(«) = r(q\Po\q) = ~<p Vr(x)e- AWx) \ 

£ £ Jr(0) = q 



(4) 



Finally, we define the partition function as a normalization factor 
which is obtained by integrating over q, 



r i f°° c r (^) = 1 
Z = Jr\ e~ m = / dq(b Vr(x)t 



-A[r(x)] 



(5) 



We are now in a position to define an equilibrium effective 
potential, which encapsulates knowledge of the physical quantum 
density into a form amenable to classical sampling. We choose the 
equilibrium effective potential, W(q) such that, 



loO?) 
W(q) 



1 



P 



log 



r(fift)= 9 



Vr(x)e 



A[r(x)] 



r(0) = q 



(6) 
(7) 



Note that this definition associates the Boltzmann factor, e^P w W, 
with the unnormalized density. Because r|o(<j) must integrate 
to unity, this allows us to easily recover the partition function 
and corresponding quantum Helmholtz free energy, A, with the 
classical integral, 



f 



dq e 



-mil) 



1' 

J—c 



dq\\ 0 (q) 



(8) 



Using Equation 7, one can exactly calculate the equilibrium 
effective potential whenever one can evaluate the path integral. 
Unfortunately that is usually numerically intractable. Thus, it is 
useful to wonder if a unique map exists between any potential 
V(q) and W(q) under the conditions of a fixed ensemble. If one 
could easily evaluate the map one could transferably adapt BO 
potentials to give physical results in classical simulations. Since 
this mapping is a functor 2 which gives an effective force-field we 
refer to the map as the "force-field functor" and denote it with 
the symbol T. A morphism depicting the structure of our proof 
is shown in Figure 1 . 



Exists (Eq. 4) 



Exists (Eq. 7) 




Unique (Section IIA) Unique (Section IIB) 

FIGURE 1 | Morphism depicting uniqueness and existence of 
mappings between the physical potential, V(q), the equilibrium 
effective potential, W(q), and the associated quantum and classical 
equilibrium densities, r\a(q) and r\o(q), respectively. This establishes the 
existence of a mapping, !F, which uniquely determines the equilibrium 
effective potential. 



throughout this paper the variable "q" refers to a position in the full coor- 
dinate space of the system (q e di 3N where N is the number of particles). 
To distinguish a position variable from a path variable we will use r(t) to 
represent a particular trajectory. 



A functor differs from a functional in that a functor maps one vector space to 
another whereas a functional maps a vector space to a scalar. In this context, 
"operator" is a more common term than "functor" but we prefer to call this 
"force-field functor theory" to evoke the connection with DFT. 
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2. UNIQUENESS AND EXISTENCE 

Our first step toward developing a theory of force-field func- 
tors is to show that the proposed mapping, .F[V(i?)] — > W(q), 
exists and is unique. This proof begins in Part A of the current 
section in which we argue that no two V(q) lead to the same 
quantum equilibrium density r|o(q), which exists by Equations 3 
and 4. To show this we take inspiration from Mermin's extension 
of the Hohenberg-Kohn theorem for finite temperatures and use 
the quantum Bogoliubov inequality to construct a proof by con- 
tradiction (Mermin, 1965). For any potential without hard-shell 
interactions, the density is always given by a Boltzmann factor 
of the potential as in Equation 6; thus, the equilibrium effective 
potential exists for any physically-relevant quantum potential. In 
Part B of the current section, we make a similar argument to prove 
that there is a one-to-one map between classical equilibrium 
density and classical potential (Chayes et al., 1984). Since the 
effective potential is chosen to be the classical potential associ- 
ated with the quantum density, these results directly imply that 
the map between physical potential and effective potential must 
be unique. 

2.1. UNIQUENESS OF QUANTUM DENSITY 

Both steps in this proof take the form of reductio ad absurdum 
arguments based on the uniqueness of an ensemble which mini- 
mizes the free energy of a canonical system. In the Appendix B we 
show that, 



we see that, 



A[p] > A [p 0 ] , p # po 
where A is the quantum Helmholtz free energy, 



A[p] = Tr 



/' ( H + i log [p] 



(9) 



(10) 



which is minimum when p is equal to the quantum equilib- 
rium density matrix po associated with the Hamiltonian, 7i = 
T+ V(q). With this in mind, suppose that there were another 
potential V(q) that led to the same density r\o(q). Denote the 
Hamiltonian, canonical density matrix and Helmholtz free energy 
associated with V(q) by H, po, and A. Since V(q) ^ V(q) and 
p Q ^ po 3 we can write 



A = Tr 



Po ( ft + p log [po. 



(11) 



Po ( ft + p log [Po. 



< Tr 

= A + Tr[p 0 V(q) - p 0 V(q)]. 
Using the definition of the quantum equilibrium particle density, 4 
r\o(q) = rr[p 0 \q}(q\], (12) 

3 That the corresponding equilibrium density matrices are not equal is obvious 
in Equation 1. 

4 Recall that q 6 9i 3N so, \q) = n?=i 



A < A + J dq [V(q) - V(q)] n 0 



(<?)• 



(13) 



But we see that this relation is still true if we interchange over- 
scored variables, 



A < A + J dq [V(q) - V(q)] n 0 



(<?)• 



This leads to the contradiction, 

A + A<A + A. 



(14) 



(15) 



and therefore only one V(q) can result in a given r|o(<j). This 
proves that V(q) uniquely determines T}rj(<j). Next, we show that 
the only potential which can reproduce the quantum density with 
classical sampling is the equilibrium effective potential. 

2.2. UNIQUENESS OF THE EFFECTIVE POTENTIAL 

Equation 7 shows the existence the equilibrium effective poten- 
tial, W(q). It remains to be shown that W(q) is the only such 
potential which will reproduce the quantum density, which is 
to say that T is completely unique. The classical Bogoliubov 
inequality states that, 

A [ifoOj)] > A [tloOz)] . ifoO?) # loOz) (16) 
where A is the classical Helmholtz free energy, 



A[x) 0 (q)] = B [%(«)] " pSM?)] 



(17) 



(q)W{q) + ^ j d<jn 0 ((j) log[n 0 ((j)] 



which is minimum when r)o(r) is equal to the classical equilib- 
rium density in the presence of W(q). For completeness, this 
result is also proved in Appendix C. With this in mind, suppose 
that there were two effective potentials, W(q) and W(q) that led 
to the same density. Then, 



A= j dq\] 0 (q)W(q) + i J dq r) 0 (q) log [f)o(q)] 



(18) 



J dq\] 0 (q)W(q) + i J dqrio((j)log [f] 0 (q)] 
A + j dqx\o(q)[W(q)-W(q)]. 



So we see that, 



A<A + y dqr\ Q (q)[W(q) - W(q)] 



(19) 



If we interchanged all over-scored quantities, we would also find 
the following, 



A < A + y dq^{q)[W{q) - W(q)] 



(20) 
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Adding these equations together leads to the result, 



on q or V(q), to write 



A + A < A + A. 



(21) 



Thus, we see that no two W(q) lead to the same r\o(q). 

Because the physical potential V(q) uniquely determines the 
quantum equilibrium density r\o(q), and the quantum equi- 
librium density uniquely determines the equilibrium effective 
potential W(q), we see that the map, T[V(q)] — >- W(q) must be 
completely unique. 

3. APPROXIMATE LINEARITY 

The results of the above section establish the possibility of revers- 
ing T by modeling pairs of V(q) and W(q) generated via the exact 
path-integral. However, the concept of J- is not useful unless we 
have good reason to suspect that T or a useful approximation 
to T will be easy to obtain and evaluate numerically. In this sec- 
tion, we analyze the approximation of T as a linear functor which 
is straightforwards to construct numerically and because of its 
separability, applicable to systems of arbitrary dimensionality. 
We begin by rewriting Equations 4 and 6, 



r(0) = q 



Vr(x)e- A[r(x)] 



(22) 



and introduce several definitions which break apart the action 
term into a kinetic part and a potential part, 



U[r(x)] = -j> dx[V(q) - V(r(x))] 

I dx m; r,( 

Jo ■ , 



T[r(x)] = exp 



1 



Zt 



f dq(f 

J -co Jr\ 



(0) = q 



Vr(x)T[r(x)]. 



(23) 
(24) 
(25) 



We now employ a notation due to Feynman and Hibbs, for the 
equilibrium average of a path functional weighted by T and nor- 
malized by Zt, "()" (Feynman and Hibbs, 1965). This allows us 
to write a concise, exact expression for W(q): 



■ Zt e 



(26) 



Jensen's inequality tells us that that, (ef) > e^ with an error on 
the order of the variance of/. This simplifies the path integral and 
introduces error that is second order at worst in the weighted path 
functional average, 



(eWLVWh _ e (U[r(x)]) 



-f,W(q) 



+ <D[(U[r(x)]) 2 - (U[r(x)] 2 )] (27) 



: Zt e 



HV(q) e Mr(T)]}_ 



(28) 



Because any potential is unique only up to a constant, we can use 
properties of logarithms to remove Zt, since it does not depend 



W(q) * V(q) - -(U[r(x)]) 
P 



(29) 



with corrections on the order oilA 2 . We also see from this that the 
equilibrium effective potential is a temperature dependent correc- 
tion to the true potential. U[r(x)] is clearly a linear functional of 
V(q) and (U[r(x)]) is clearly a linear functor of U[r(x)], 

I rr(m = q 

<W[r(T)]> = — 4 Vr(x)T[r(x)]U[r(x)]. (30) 

ZT Jr(0) = q 

In the multi-dimensional case, the path integral couples all 3N 
modes of q, making the exact T a very complicated object which 
embeds all-orders of quantum many body effects between these 
modes. However, our analysis suggests a linear approximation 
which conserves the locality of the original potential. With this 
approximation we can separate the integral in Equation 30 into 
each individual interaction order of the potential and see that the 
path integral does not multiply these terms; the pairwise interac- 
tions remain pairwise, the three-mode interactions are mapped 
by T onto three-mode interactions, etc. 

Approximate separability of this mapping is one of the key dif- 
ferences between our method and approaches such as Feynman- 
Kleinert, which introduces higher ordered many-body terms into 
the effective potential, or RPMD, which avoids the issue at the 
cost of introducing ancilla particles. Our T can be imagined as 
a Gaussian smearing of V(q) to first approximation. It is rea- 
sonable to suspect that the non-separable many body couplings 
would be blurred to a high order such that the many-body expan- 
sion of the equilibrium effective potential might terminate faster 
than the many-body expansion of the uncorrected physical poten- 
tial. This agrees with the empirical observation that tunneling 
effects stabilize pairwise interactions more than higher-ordered 
interactions. 

4. NUMERICAL TESTS 

It is far from obvious that a transferable map between V(q) and 
W(q) can be practically obtained and usefully accurate. Instead of 
calling upon the most sophisticated procedures we can implement 
to solve the problem, we take the simplest approach to devel- 
oping and testing our approximation to T so that our results 
are designed to be a worst-case, upper-bound on the error that 
leaves room for optimism. Approaches such as machine learn- 
ing could be employed in future work (Snyder et al, 2012). We 
approximate J 7 as a linear map (a matrix) acting on our potentials 
vectorized into coefficients of Legendre polynomials. The entries 
of this matrix are determined by simple least-squares on a ran- 
domly generated training set of 1000 one-dimensional potentials 
and their corresponding effective potentials chosen by randomly 
choosing Legendre coefficients with the only constraint being that 
the classical densities vanish at their boundaries. 

Effective potentials were calculated using Equation 7 with den- 
sities obtained from the efficient real-space discrete variable rep- 
resentation (DVR) of the path integral (Thirumalai et al., 1983). 
We examine how this T performs on instances of other random 
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FIGURE 2 | Top: plot of the percent error in potential energy of a classical 
simulation with the classical potential (blue) and F generated distribution 
(red) against Q. Bottom: plot of mean integrated squared error (MISE) from 
the exact quantum density for classical density (blue) and !F generated 
density (red) against Q. Each point is the mean of these errors on 1000 
random potentials with 50 basis functions. 



potentials not included within its training set and then apply it 
to the Silvera-Goldman pair potential for liquid para-hydrogen 
(Silvera and Goldman, 1978; Nakayama and Makri, 2003; Hone 
and Voth, 2004; Poulsen et al., 2004; Miller and Manolopoulos, 
2005). Using the resulting effective potential, a classical Monte 
Carlo simulation was performed to give us radial pair distribu- 
tion functions in agreement with results from PIMC at a fraction 
of the computational cost. 

4.1. OBTAINING THE LINEAR FUNCTOR 

In order to obtain the simplest possible T we model the linear 
transformation as a matrix. This requires that we treat the phys- 
ical potential and effective potential as vectors in some basis of 
real-valued functions. Because force-fields are often chosen for 
the speed with which they can be evaluated, it seems natural 
to use a polynomial basis. Legendre polynomials evaluated on a 
fixed domain of [—1, 1] were chosen for their orthogonality and 
historical usefulness in fitting potentials. 

Consider the short time Trotterization of the path-integral, 
which we use to generate exact quantum densities for our test sets 
(Thirumalai et al., 1983). The short-time propagator effectively 
acts as a Gaussian which blurs out the density with a variance 
that depends exactly on the inverse of the square root of the 
the mass times the temperature. This factor which determines 
the "quantumness" of the system is proportional to the thermal 
de Broglie wavelength, A = h^/2nfi/m (Yonetani and Kinugawa, 
2003; Georgescu et al., 2013). Because we wish to calculate the 
deformation of a potential vector evaluated on a fixed domain, 
the parameter which characterizes our map must depend on the 
ratio between the thermal de Broglie wavelength and the potential 
length-scale, Q = A/L where L is the potential length-scale. 

In order to obtain a linear functor capable of transform- 
ing a one-dimensional potential at fixed Q into another one- 
dimensional potential at fixed Q we randomly generated pairs of 
potential vectors and their corresponding effective potential vec- 
tors. These vectors were in a Legendre polynomial basis of order 
B and the vector elements of the classical potential (i.e., basis 
coefficients) were drawn from a flat distribution between — 10/(5 
and 10/(5. The corresponding effective potential vectors were cal- 
culated by evaluating the classical potential vectors as Legendre 
polynomials on the fixed domain, passing the scalar potential 
and Q to the aforementioned DVR routine which yields a scalar 
quantum density, and finally fitting the negative logarithm of that 
density divided by (5 to a vector of Legendre polynomials in accor- 
dance with Equation 7. Having done this, the goal is to find a 
matrix T € B x B such that, TV ~ W. We chose to perform a 
Levenberg-Marquart L2 optimization to determine the elements 
of this matrix (Levenberg, 1944). Our residual was defined as 
the concatenation of the difference vectors, TV; — Wj for all 
N physical potential / effective potential pairs in the randomly 
generated set. 

4.2. PERFORMANCE ANALYSIS 

The linear approximation to T appears to work quite well for 
even fairly large values of Q. As we can see in Figure 2, the errors 
on an independent test set from the linear T generated W(q) 
are minimal and significantly better than the classical predictions, 



especially in strongly quantum regimes. Even the deviation from 
the exact answer is improved relative to simulations which employ 
the uncorrected physical potential. For both simulations the error 
goes to zero as Q goes to zero — a consequence of classical cor- 
respondence. As one might expect as Q is increased, predictions 
given by both classical and T generated distributions deviate 
more significantiy from the exact answer. In the W(q) simulations 
these errors are entirely due to the linearity of T. Another view 
of the the performance of the linear functor is given in Figure 3. 
When temperature and length are fixed, mass is a reasonable pre- 
dictor of the performance of both W(q) and V(q) simulations. 
For low masses, the classical treatment often misses the quan- 
tum free energy by as much as a kcal/mol (chemical accuracy). 
Having characterized the error of assuming linearity we turn to 
separability. 

Figure 4 shows the effective potentials obtained from apply- 
ing our linear T, trained at 14K and 25L with sets of 1000 
potentials, to the Silvera-Goldman potential, which is perhaps 
the most common potential used to simulate liquid hydrogen 
with path integrals (Silvera and Goldman, 1978; Nakayama and 
Makri, 2003; Hone and Voth, 2004; Poulsen et al, 2004; Miller 
and Manolopoulos, 2005). We then performed a classical Monte 
Carlo simulation on the potential mapped at 25K and the poten- 
tial mapped at 14K, using 150 molecules in a cubic cell with 
periodic boundary conditions and one million steps. Cell size 
was fixed by densities from the literature (Nakayama and Makri, 
2003). 
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FIGURE 3 1 Left: correlation of classical (blue-green) and F generated 
(red-yellow) free energy with exact free energy. Dotted lines enclose 
the chemically accurate region of within one kcal/mol. In more than 
97% of instances, our map is more accurate than the classical 
treatment. Right: correlation of classical and F generated potential 



energy with exact potential energy. Color brightness indicates the 
mass used in setting the Q value at 25K. As mass increases, 
classical simulations better approximate the energy. Data consists of 
1000 cross-validating potentials at each of the six masses shown on 
the colorbar. 



3.0 




Radius (Bohr) 

FIGURE 5 | The top box shows radial pair distribution functions at 
25K and the bottom box shows radial pair distribution functions at 
14K. The blue (dotted-dashed) curve is for the classical liquid without 
correction for quantum effects. The green curve (solid) shows the result 
of classical Monte Carlo sampling on the effective potential obtained 
with our linear The red curve (dashed) shows PIMC results 
(Nakayama and Makri, 2003). Even this simple :F is a major 
improvement over the classical potential. 




I I I i I I I 

4 5 6 7 8 9 10 11 12 

Radius (Bohr) 



FIGURE 4 | The dashed black line above shows the classical 
Silvera-Goldman potential in the region of interest for our problem. 

The red line is the effective potential obtained with our linear J 1 at 25K and 
the blue line is at 14K. 



The resulting radial distribution functions, g(r) are shown 
in Figure 5. The differences between the W(q) generated g(r) 
and the PIMC results are presumably due to the assumption 
of separability. Slight over- structuring of g(r) at the first shell 
results from neglect of the 3-body components of the exact W(q). 
Remarkably, this over-structuring appears to decrease with tem- 
perature, lending credence to the idea that many-body effects 
in W(q) are largely blurred-out by the smearing which the low 
orders of T perform on the potential. At both temperatures the 
errors of these approximations are quite reasonable and although 
the classical system undergoes a non-physical transition to a solid 
between 25 and 14K, the model of the present work remains 
correct. 

5. CONCLUSION 

We have shown that for each physical potential, there is a unique 
effective potential which reproduces the quantum density and 
free energy when sampled with classical statistics. Other prop- 
erties of a classical simulation of the effective Hamiltonian are 



not designed to approximate reality by the mapping, but the 
effective potential maybe advantageous to the status quo: classical 
simulation on a Born-Oppenheimer surface. In this paper we 
have shown that the implied mapping between the physical and 
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effective potential, T , can be made concrete to a useful degree 
of accuracy. A simple linear model for T improves on the 
physical potential systematically over a broad range conditions. 
Even under the assumption of separability and without any 
exponential functions in our training set, our model for T 
adequately describes the density of a popular para-hydrogen 
model at exactly the cost of the corresponding classical simu- 
lation. Non-linear models for T and expressions which do not 
assume complete separability are likely to improve on these results 
and produce even more accurate transferable recipes for digest- 
ing Born-Oppenheimer potentials. In particular, we imagine the 
development of simple functors which could be applied to Born- 
Oppenheimer surfaces so that classical sampling will immediately 
give improved results. Ultimately, we believe that force-field func- 
tors can provide a scalable methodology for including quantum 
propagation effects in systems that are intractable for exact meth- 
ods, such as protein force-fields. 
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A1 . APPENDIX We see that p x = po if A- = 0 and p x = p if X = 1. Accordingly, 
A1.1. QUANTUM DENSITIES FROM CLASSICAL SAMPLING 

In practice, path integral expressions are analytically intractable , r , , r , f 1 3 , . , , 

. / ~ . ■ • r A[p]-A po = / — A[px]dX (A8) 

except in a tew cases. Feynman proposed to simplify Equation 5 J 0 dX 

by changing from an integral over all closed paths that start and 

end at point q to an integral over all closed paths that have an by the fundamental theorem of calculus. To evaluate the deriva- 
average value equal to the path centroid r, tive we use, 



1 /-P* 

- / dtr(t). (Al) 
Jo 



A[p x ] = (A9) 



Tr 



,,, [ H + XA + ~log[p k ] 



XTr[Ap 



So that we only integrate over each closed path once, we must 
change our expression for the partition function to only calculate 

paths that match the centroid, The ^ rst trace is stationary for variations of p x about the corre- 

sponding density matrix. Thus, we only need to differentiate the 
second trace, 



Z= j dq<j)Vr(T)h[q-r]e 



A[r(x)] 



(A2) 



3 

— A [ml = — XTr 
Vr(z)e- A[r(x)] . dX LWI 



3 

A — p x 
dX 



(A10) 



While the partition functions given by Equation 5 and We evaluate m p x using the operator identity, 
Equation A2 are exactly equal, the two expressions are associated 

with subtly different scalar density functions. Equation 5 is asso- 3 p -p(H+ X A) _ g -$(H+lA) f ^ ^ (All) 

Jo 



ciated with the true equilibrium density in Equation 4 and A2 is dX 
associated with the path centroid density, 



1 



0 

A, (P') = e P'(«+^) Ae -P'(W+^A) (A12) 



ipr(T)5[q-f]r^ [r(,)l . (A3) 3 /-P , r . ,. 

J -p x = -J^dVp l [A x (V)-(A) l ], (A13) 

The Dirac delta function in this equation enforces the require- 
ment that integrating the Boltzmann factor associated with this where 

density over the path centroid, r, will result in exactly the path (X) x = Tr[p^X] . (A14) 

integral expression for the quantum partition function (Jang and 

Voth, 1999). The centroid density plays a prominent role in CMD Therefore, 

and Feynman-Kleinert methods but does not apply to force-field 

fUnCt0rthe ° ry - ^A[p x] =xfdV( { AA x (V)) x - { A)l). (A15) 

A1.2. PROOF OF QUANTUM B0G0LIUB0V INEQUALITY ° 

The quantum Bogoliubov inequality is proved for the grand By cyclically permuting operators within the trace, one can verify 
canonical ensemble in the Appendix of (Mermin, 1965). We adapt ma t 
this proof for the canonical ensemble, in the interest of com- 
pleteness, to show that for all positive definite p with unit trace, ^ x (ft)) X = (A) x Vf5', (A16) 

A[p] > A [po] , p £ po (A4) 

if A is the quantum Helmholtz free energy of the canonical 
ensemble, 

A[p] = Tr p(n+^log[p]j , (A5) 



(AAx(p')) x = ^A x Qp'j AxQp'jj. (A17) 
With these identities, we can rewrite Equation A15, 



which is minimum only when p is equal to the quantum equi- —A [p x ] = (A18) 

librium density matrix po associated with the Hamiltonian, TL = dX 
T + V(q). To start we define, ,,p 



P X =e 



H( Al G e ' 



X / dp' A X -p' - (A) x ) ( A x I -P' I — (A) 



- m -xA) /Tr ^ e - m+ xA^ (A6) 

A = l0g[p] -H. (A7) " r -". f T ~V r " r- --— r 

P &L ^ J v ' minimum ot the tree energy must occur when p x = po. 



w ^ ere This integral is non-negative and can be zero only if A is a mul- 

1 tiple of the unit operator, i.e., if po = p. This proves that the 
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A1.3. PROOF OF CLASSICAL B0G0LIUB0V INEQUALITY 

If r|o (q) is the equilibrium density for a classical canonical ensem- 
ble and rfo(q) is a different density, Gibbs' classical Bogoliubov 
inequality states that, 

A [r\ 0 (q)] > A [r] 0 (q)] , r\ 0 (q) # r) 0 (<?) (A19) 

where A is the classical Helmholtz free energy, 



A[t]o(?)] =E[r) 0 (q)]- -S [tlo(q)] 



(A20) 



/d, 



Tio((j)W((j) + i y dqt) 0 (q) \og[^(q)] . 



To see that this is the case we start by writing, 

i j dqTf 0 (<j)log[ifo((j)] > dqtf 0 (q)log[tio(q)]. (A21) 
The difference between the right and left sides of this equation is, 

^ J &q (ifoC?) log [rfoO?)] ~ r\o(q) log[rio(q)]) (A22) 

i f , - , ri~o(q)" 

= 7; / d<JTlo(q)log — — . 

Because log [x] > 1 — ^ and we know that the densities are 
normalized, 

\ f dqif 0 (<?)log -If d <l[no(q) - r\o(q)] =0. 

(A23) 

We can simplify this further to, 



i log [rio(q)]\ > /i log [rfo(<?)] 



(A24) 



We know that, 



TloOj) 

nbOz) 



e -P(?) 
z 



where £(q) and Z correspond to the energy and partition function 
associated with r\o(q). Thus, 



log 



-E(q) 



z 



-log [Z] 



(.log 



z 



(A27) 



-£(<?) - - log [Z] 
P 



We may safely assume that (E(q)) = so using the definition 

of the Helmholtz free energy, A = — ^ log[Z], 



A [rfo(q)] > A [tio((?)] , rfoOj) 7^ t)o(<?)- 



(A29) 



A1.4. APPLYING LINEAR FUNCTOR TO SILVERA-GOLDMAN 

The matrix which was ultimately used to transform the Silvera- 
Goldman potential was obtained by fitting 1000 random poten- 
tials with B = 50 basis functions in the appropriate Q regime. The 
Silvera-Goldman potential has the form, 



V(r) = exp[ct — 8r — yr 2 ] 

/ C 6 C 8 C10 \ C 9 



(A30) 



(r) 



where 



/ c (r) 



e-^A-i)^ if r < rc 
1, otherwise. 



(A31) 



Parameters for the Silvera-Goldman potential are provided in 
Table Al (Silvera and Goldman, 1978). 

Exponential functions cannot be easily represented in a poly- 
nomial basis and the Silvera-Goldman potential diverges expo- 
nential as r approaches zero. Accordingly, we fit the potential only 
in the physically relevant region of r > 4 Bohr. We matched the 
slope of the potential at r = 4 Bohr and extend the potential as 
a straight line in the region 0 < r < 4 Bohr. We choose to fit 
the potential out to r = 24 Bohr but imposed a standard cut- 
off after the fact at r = 20 Bohr as the potential is clearly flat 
by this point. We simulated para-hydrogen at 14K and 25K. At 
25K, the thermal de Broglie wavelength is 4.6 Bohr; thus, a cut- 
off distance of 20 Bohr gives Q = 0.23. At 14K, the thermal de 
Broglie wavelength is 6.2 Bohr and Q = 0.31. Based on statis- 
tics collected from 10,000 random potentials generated with these 
Q values, in both of these regimes, the classical free energy is 
more accurate than the T predicted free energy less than 1% of 
the time. 



(A25) 

Table A1 | Parameters of the Silvera-Goldman potential (Silvera and 
(A26) Goldman, 1978). 



Parameter 


Value (atomic units) 


a 


1.713 


8 


1.5671 


Y 


0.00993 


Ce 


12.14 


C B 


215.2 


c 9 


143.1 


C10 


4813.9 


re 


8.321 
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